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We study an open-boundary version of the on-off zero-range process introduced in Hirschberg et 
al. [Phys. Rev. Lett. 103, 090602 (2009)]. This model includes temporal correlations which can 
promote the condensation of particles, a situation observed in real-world dynamics. We derive the 
exact solution for the steady state of the one-site system, as well as a mean-field approximation for 
larger one-dimensional lattices, and also explore the large deviation properties of the particle current. 
Analytical and numerical calculations show that, although the particle distribution is well described 
by an effective Markovian solution, the probability of rare currents differs from the memoryless case. 

In particular, we find evidence for a memory-induced dynamical phase transition. 
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I. INTRODUCTION 

The first step in the study of a complex system is often 
to focus on its typical behaviour. Employing the tools of 
statistical mechanics, for example, we can study the typ¬ 
ical properties of a macroscopic system. However, there 
are situations in which the behaviour of interest is not 
typical, but rather atypical. For example, the transport 
of energy, particles or vehicles could be enhanced by ex¬ 
ceptional coherent configurations, or occasionally delayed 
when an instantaneous situation similar to condensation 
occurs [Tj. As another example, in communication net¬ 
works it is very important to predict how likely it is to 
have interruptions or packet loss [2]. 

Moreover, rare events help to shed light on the foun¬ 
dations of non-equilibrium statistical mechanics, just as 
they play an important role in defining the thermody¬ 
namic potentials in equilibrium statistical mechanics 0- 
[5j. While a general framework for the characterisation 
of systems far from thermal equilibrium is at a primitive 
stage, large deviation theory plays a central role 0®. 

As a comprehensive theory of non-equilibrium phe¬ 
nomena does not exist, the analytical study of toy models 
is an effective way to build it up. The majority of the 
interacting-particle models in the literature are Marko¬ 
vian, i.e. memoryless. Such an approximation simplifies 
the theoretical treatment, but can exclude some proper¬ 
ties of physical phenomena. The effects of memory on 
such models have thus prompted recent curiosity [2HlQ|. 
We enter into this context by studying a driven-diffusive 
system which is referred to as the on-off zero-range pro¬ 
cess (on-off ZRP) and focusing on its non-equilibrium 
aspects. This model is an open-boundary version of the 
non-Markovian ZRP defined in Hirschberg et al. CO ® 
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and allows analytical progress. 

Non-equilibrium stationary states (NESSs) are charac¬ 
terised by the presence of finite currents, which measure 
the violation of detailed balance for opposing transitions 
between two configurations dHE]. Such currents fluc¬ 
tuate in time and the functions that determine the proba¬ 
bility of deviation from their typical values have the same 
mathematics as the thermodynamical functionals defined 
in equilibrium statistical mechanics. In this spirit, we are 
interested in the particle current for our model, i.e., in 
the transition events corresponding to particle hops. The 
study of its rare fluctuations reveals effects of the time- 
correlation hidden in the stationary state. 

The paper is organised as follows. In Sec. [TTJ we define 
the model and in Sec:. |1 11 A] we derive its stationary state 
in the single-site system. In Sec. |IIIB we present a mean- 
field treatment of the dynamics on a chain topology and 
test the validity of this approximation against standard 
Monte Carlo simulations. In Sec. |IV[ we explore the cur¬ 
rent fluctuations, focusing, in a one-site system, on the 
difference between the small fluctuation regi me, obt ained 
by analytic continuation of the NESS ( Sec. IV B 1 1, and 
the extreme fluctuation regimes (Sec. IVB3I, and de - 
riving the phase boundaries between them (Sec. IVB2|. 
The analytical results are tested against an advanced nu¬ 
merical method which has been developed to evaluate 
large deviation functions directly na. We summarize 
the results in Sec. [V] 


II. MODEL 

The ZRP is a model of interacting particles on a dis¬ 
crete lattice, which we take here to be a one-dimensional 
chain. Each lattice site can contain an arbitrary positive 
number of particles. The evolution proceeds in continu¬ 
ous time, i.e., transitions occur after a waiting time which 
is an exponentially distributed random variable. Specif¬ 
ically, in the standard ZRP, a particle can hop to one 
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of the neighbouring sites with rate proportional to p n , 
which depends only on the occupation number n of the 
departure site. Obviously, the departure rate from an 
empty site is given by /To = 0. The functional form of p, n 
encodes the interaction between particles, which occurs 
only on the departure site, hence the epithet zero-range. 
The special case /z ra = an, where a > 0 is a constant, 
corresponds to free particles since in this case each parti¬ 
cle leaves the site independently from the others. Other 
choices of /r n correspond to attractive or repulsive inter¬ 
particle interaction if the n-dependence is sublinear or 
superlinear, respectively. 

Models with zero-range interactions have proven to dis¬ 
play complex collective behaviour whilst allowing analyt¬ 
ical treatment m- In particular, the ZRP is well suited 
for theoretical analysis because the stationary probabil¬ 
ity distribution of a given configuration factorises and 
can be calculated exactly jlfij . It is worth mentioning 
that certain choices of ji n lead to condensation, i.e., the 
accumulation of a macroscopic fraction of the total num¬ 
ber of particles on a single site. Condensation transitions 
far from equilibrium have been studied in physics m , as 
well as in economics mm, biology [2D], network sci¬ 
ence num, and queueing theory [23]. Toy models, such 
as the ZRP, provide a theoretical foundation for under¬ 
standing condensation in these systems. Exact results 
from the ZRP have also been used in models of vehicu¬ 
lar traffic [U [23] , reptation in polymer physics (2H3, and 
transport and coalescence in granular systems | 26 j . 

A further step towards a deeper comprehension of real- 
world phenomena may be achieved by studying stochastic 
systems with time correlations. A modified zero-range 
process with non-Markovian dynamics has been intro¬ 
duced in Hirschberg et al. Hi- The crucial new in¬ 
gredient is that each site has an additional clock/phase 
variable r and the particles cannot leave the site when the 
clock is set to zero, which corresponds to the OFF phase. 
The clock ticks and turns on with rate c and turns off 
with each particle arrival. This mechanism favours the 
accumulation of particles on a site. According to the 
zero-range dynamics, the particles interact only on site, 
but now have a different departure rate /x raiT . The ad¬ 
ditional variable r takes into account events in particle 
configuration space that happened in the past and there¬ 
fore introduces temporal correlations. This model has 
sparked interest as it displays, under certain conditions, 
a condensate with slow drift [3 IS- Systems with dis¬ 
tinct on and off phases are also of interest as models 
for intra-cellular ion-channels EHIMi and for data traffic 
streams [2]5], as well as providing examples of stochastic 
processes with non-convex rate functions m- 

To the best of our knowledge, the ZRP with on-off 
dynamics has been studied only on ring topology, i.e., 
with periodic boundary conditions. In this paper we in¬ 
vestigate the open-boundary version of the model, thus 
extending the work of Hirschberg et al. We imple¬ 

ment the same dynamics on an open chain with arbitrary 
hopping rates and boundary parameters, see Fig. [T] Par- 
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FIG. 1. (Color online) Non-Markovian ZRP on a one¬ 
dimensional lattice with open boundaries. Each site has a 
hidden variable r, whose values are represented by the po¬ 
sitions of a gear, which controls the departure rate. When 
r assumes value zero no departure is possible and the corre¬ 
sponding state is referred to as OFF. This lock-down occurs 
in conjunction with the arrival of a particle. 


tides are added and removed through the boundaries. 
On the leftmost lattice site (site 1), particles are injected 
with rate a and they are removed with rate 7/t n ,T which 
is non-zero only when the phase of site 1 is different 
from r = 0. Similarly, on the rightmost site (site L) 
particles are removed with rate /3/T niT , according to the 
phase of site L, and are injected with rate S. This situ¬ 
ation corresponds to a bulk system in contact with two 
different reservoirs. In the bulk, particles jump to the 
left (right) with rate qn U}T (p/r njT ), which is again non¬ 
zero only when the phase of the departure site is not 
r = 0. The dynamics is sensitive to the specific rate 
values and we consider choices that induce a rightwards 
driving along the chain. In particular, it is worth making 
the distinction between the totally asymmetric (TA) and 
the partially asymmetric (PA) processes. 

Hereafter, we will consider explicitly two forms for in¬ 
teraction factor pLji,r , i.e., the case where /r„, r is constant 
with respect to n > 0, which corresponds to an on-site at¬ 
tractive interaction between particles, and the case where 
Hn }T is linear in n, which corresponds to no direct inter¬ 
action between particles (excluding residual correlations 
due to the blockade mechanism). 

The stationary distribution of the standard ZRP on an 
open chain has been extensively studied in Levine et al. 
ED- In this case the particles are distributed along the 
system according to a product-form structure that im¬ 
plies no correlations between sites. In contrast, the on-off 
ZRP can generate more complex patterns, as shown in 
Fig- a for three sets of parameters. The clock-tick rate c 
plays a major role in these patterns, the lower its value 
is, the more important the correlations are. Increasing 
the value of c, the system eventually becomes spatially 
uncorrelated. In the next two sections we study in de- 
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tail how the introduction of time correlations affects the 
stationary state and the current fluctuations. 


III. STATIONARY STATE 
A. Exact results for one-site system 

As mentioned above, a notable property of the stan¬ 
dard ZRP is that in the stationary state the probabil¬ 
ity P*({m}) of finding the system in a state {n;} = 
(ni, ri 2 , ■ . ■ ul) with ni particles on site l, is given by a 
simple factorised form 


(a) 



L 

P*({ni}) = l[Pnn l ), (1) 

1 = 1 

where P*(n/) is the probability of finding the site l with 
ni particles. The one-site marginals are determined by 

n n 

W = yII^ ( 2 ) 

where zi is a site-dependent fugacity (which is a func¬ 
tion of the hopping rates) and the grand-canonical par¬ 
tition function Z/ = Zl " 11"= l H i 1 ensures normal¬ 

isation [3T]. It is worth noting that, for certain choices 
of n n , it is not possible for the sum in Zi to converge 
for all zi. The divergence of Zi corresponds to the accu¬ 
mulation of particles on the site l and we refer to it as 
congestion. Indeed, the infinite accumulation on one or 
more sites in an open system can be thought of as a kind 
of condensation phenomenon H3GH- In the following, 
we will also use the “condensation” terminology even for 
the single-site case. 

Our preliminary simulations in Fig. [2] suggest that we 
cannot rely on a factorised steady state for the non- 
Markovian model introduced in Sec. [TTJ However, for the 
single site system, an exact solution is straightforward. 
The state is defined by two variables: the number of par¬ 
ticles in the box n and a “clock” or “phase” variable r. 
We focus on TA dynamics and consider a box which re¬ 
ceives particles with rate a and ejects particles with rate 
(3Hn,T , where /r n , T is a function of the box state. The 
departure event is possible only when r ^ 0. Also, the 
dynamics includes the advance of the clock with rate c, 
and the reset to r = 0 when a particle arrives. If one de¬ 
fines P(—l,T;t) = P(n, — l;t) = 0, the following Master 
equation is valid for r > 0 and n > 0: 

dP(yn 'J' ^ = cP(n, t - 1; t) + Pn n+l T P(n + 1, r; t) 
at 

+ S T , o y] aP(n — 1, t'; t) — (c + /3 /x„ jT + a)P(n, r; t), 

t’> 0 

( 3 ) 

where P(n,T\t) denotes the probability of finding the 
system with n particles and phase r at time t and S t ,q is 


(b) 




FIG. 2. (Color online) Monte Carlo time evolution of the 
occupation profile of the on-off ZRP on a one-dimensional 
lattice, (a) Rates fi n ,T = n if r > 0, /r n ,r = 0 otherwise, and 
(ot,/3,y,6,p,q,c) = (0.1,0.2,0,0,1,0,0.05). Only correlations 
due to the blockade mechanism are present. The particles 
organise in travelling clusters. Their speed is mainly governed 
by c. (b) Same parameters as the former case, except (/3,p) = 
(10 4 ,10 4 ). The particles jump almost simultaneously to the 
next site as soon as the blockade is removed. Each particle 
cluster tends to occupy a single site. The drift proceeds with 
a rate ~ c. (c) Same parameters as (a), except p, ntT = 1 if 
r > 0 and n > 0, p n ,r = 0 otherwise, and c = 0.15. As a 
result of the attractive inter-particle interaction, the clusters 
with more particles travel slower than the less populated ones. 
This mechanism enhances congestion. 
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a Kronecker delta. The first term on the right-hand side 
of ([3]) corresponds to a clock tick, the second term to the 
departure of a particle, the third term to the arrival of 
a particle and the fourth term to the respective escape 
events from the state (n,r). 

As in Hirschberg et al. we choose to simplify the 

dependence of the jump rate on r to /r niT = fi n when 
r > 0. Hereafter we specialise to this case, except when 
we explicitly refer to a general form for /^„ >T . In this 
simplified case, it is more convenient to write the Master 
equation © in terms of P(n, ON;*) = E r >o P ( n ’ T > 
and P(n , OFF; *) = P(n , 0; *) 


dP(n, ON; t) 
dt 

dP(0, ON; t) 
dt 

dP(n, OFF; t) 
dt 

g?P(0, OFF; t) 
dt 


=cP(n , OFF; t) + (3ii n+1 P{n + 1, ON; t ) 


- Wun + a)P{n, ON; *), 


(4) 


=cP(0, OFF; t) + PmP(l, ON; t) 

— aP(0, ON; t), (5) 

=aP(n — 1, ON, t) + aP(n - 1, OFF, t) 

- (c + a)P(n, OFF;*), (6) 

= - (c + a)P(0,OFF;*). (7) 


By equating the left-hand sides of Eqs. @ 0 to zero, 
we find that the stationary distribution is given by 


P *( n ) =^~W_ W ch 


P*(n, OFF) = 


P/J"n 


a + c + /3 11 , 


-P>), 


P*(n, ON) = < a + C) P*(n), 

OL + C + p [l n 


( 8 ) 

(9) 

( 10 ) 


where z = a//3, w C}i = Hi{a + c)/{a + c + j3/ii), 
Zc = Er=o^ n nr=i w M> and = P* in, ON) + 

P*(n, OFF) by construction. We recognise the same sta¬ 
tionary state ([ 2 ]) as the standard ZRP, with an effective 
departure rate u> CjTl = /i„P(ON|n), where P(ON|n) = 
{a + c)/(a + c + f3fi n ) is the conditional probability of 
finding the site in the ON state, given that there are n 
particles. For c —> 00, the effective jump rate converges 
to the microscopic rate, i.e., w c ^ n —> ii n . The station¬ 
ary probability distribution of the occupation number is 
checked against Monte Carlo simulations in Fig. [3] for 
both constant and linear departure rates. Its tail is longer 
than the corresponding Markovian model (c — ► 00). The 
derivation of (|8|-(10) is reported in Appendix |A| 

The normalisation condition P*(n ) = 1 

on the probability distribution ([8]) requires 
limn^oo a/(/3w c P < 1. For /r n = fi, the effective 
departure rate is referred to as w c , and this stationarity 
condition is simply a/(/3w c ) < 1. It implies that values 
of c smaller than the threshold C\ := a 2 /(/3/r — a) exclude 
any stationary state and produce a congested phase. 
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FIG. 3. (Color online) Occupation probability distribution 
of the one-site system for constant ( fi n = /i,n > 0) and lin¬ 
ear (/r n = n) microscopic departure rate. The arrival and 
departure rate are a = 0.1 and /3 = 0.2 respectively. 


The onset of congestion in a larger system with constant 
departure rate is explored in Sec. 01 B| using a mean-field 
approach. For unbounded microscopic departure rates, 
i.e., limji^oo fi n = oo, the effective interaction is still 
bounded, since linin^oo w c ^ n = (a + c)//3. However, as 
long as c > 0, the normalisation condition is always 
ensured. Obviously, the linear departure rate case falls 
into this category. 

We now outline the quantum Hamiltonian represen¬ 
tation [32] of the Master equation @-0, which will 
turn out to be convenient for the study of fluctuations 
in Sec. IV The epithet “quantum” has become standard 
in the literature on interacting-particle systems, along 
with the warnings that underline that the generator of a 
stochastic process is in general non-Hermitian, contrary 
to the operators in quantum mechanics. In this approach 
one works in the joint occupation and phase vector space, 
defining a probability basis vector | n,r) = |n) G> |r) 
representing a configuration with n particles and phase 
r. A probability vector |P(*)) = En,T P ( n > r > ^)l n > r ) 
obeys the normalisation condition (1|P(*)} = 1 where 
(1| = (1,1,...). The Master equation then reads as 


j t \Pit)) = -H\P(t)), (11) 

where the operator H is a single-site Hamiltonian. Our 
convention is to use the basis kets (1,0) T and (0,1) T for 
the states |OFF) and |ON), respectively. A configuration 
with n particles is represented by a basis ket with the n- 
th component equal to 1 and the remaining components 
equal to zero. Consequently, the Hamiltonian is written 
as 


H = -c(a^ i -gT 1 )-a{a^jT 1 -t)-P{a Ni dT 1 -d Nl d Tl ), 

( 12 ) 
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with 


CLrp - 1 I 


<77i — 1 1 


( 1 0 )’ fT ' = 1 


1 0 
0 0 


, dr i=l' 


1 1 
0 0 

0 0 
0 1 


(13) 


a N 1 ~ 


l N 1 


rfjVi = 


/ 0 0 0 ... \ 
10 0 
0 1 0 

v ; 

/0 m 0 0 

0 0 H 2 0 

0 0 0 H 3 

0 0 0 0 


V: 

/O 0 0 0 

0 0 0 

0 0 H2 0 

0 0 0 p 3 

Vi 






> 1 . 


(14) 


(15) 


(16) 


We use the convention that a ladder operator with sub¬ 
script Ni or Tj acts non-trivially only on the occupation 
or phase subspace respectively. The additional subscripts 
1 underline that this Hamiltonian generates the dynam¬ 
ics for the single-site case. The operator H has a block 
tridiagonal structure which occurs in general in stochas¬ 
tic generators of processes with two variables, n and r 
in this case. The blocks correspond to changes in the 
first variable, while the entries inside the blocks corre¬ 
spond to changes of the second one. All the variables 
can change by at most 1. Such processes belong to the 
class of quasi-birth-death processes and are simple cases 
of queues with Markovian arrival and general departure 
law [33) 131 ] ■ We mention also that the results in this sec¬ 
tion can be adapted to the more general PA case with the 
replacement a —> a + 6 and /? —> (3 + 7 . Specifically, the 
quantum Hamiltonian for PA dynamics on a single-site 
is 


H = -c{a^ i -g Tl )-<x{.aN i fT 1 -t)-P(a Ni dT 1 -d Nl d Tl ) 
- - d Nl d Tl ) - S(ajfJ Tl - !)■ (17) 


B. Mean-field solution for the L-site system 

In the case considered so far, particles arrive on the 
site from the boundaries according to a Poisson process. 
The many-site system is rather more complicated than 
this. In fact, each site receives particles according to a 
more general point process, which alternates time inter¬ 
vals with no events (corresponding to the OFF phases of 
the neighbour sites) and periods with arrivals. Moreover, 


the exact statistics of the phase switching is not a priori 
known. 

In this subsection, we derive an approximate solution 
for the stationary state of the on-off ZRP on an open 
chain. The approximation consists in decoupling the 
equations which describe the dynamics for each site, re¬ 
placing the point process that governs the arrival on each 
site by a Poisson process with an effective characteristic 
rate. The decoupling of the equations allows us to us e 
the results obtained for the one-site system (Sec. |III A ). 

Let us first consider the general model described in 
Sec.[TT| where the departure rates /x„ iT retain a non trivial 
dependence on both n and r. We assume a product mea¬ 
sure Y\i=i P*( n h T l) l° r th e joint probability P*({rq, t;}) 
that the system is in a steady state with the generic site 
l in configuration (rq,T;). Imposing this solution in the 
stationarity condition of the L-site Master equation, we 
get 


cP*(ni,Ti - 1 ) +p'^2zi- 1 P l *(m - 1 , t)6 TIi o 

T 

+q' 52 z i+ iP i*( n i- 1 , T ) S T l , o +(p+q)vn l + i ,T l P*(ni+ 1 , T i ) 

r 

- [pzi -1 + qzi +1 + (p + q)Pm,n + c \P*( n i> T l) = 

(18) 

for the generic bulk site l, 1 < l < L. We use the sym¬ 
bol zi , already adopted in Sec. |III A| for the fugacity, to 
denote the ensemble average of the departure rate, since 
zi = En^^nPrM. The use of an average in¬ 
teraction term justifies the appellation mean-field. Simi¬ 
larly, for the leftmost and rightmost sites we get, respec¬ 
tively, 


cPl (n lt n - l) + a y Lj *(m - 1 , t)5 Ti ,0 

T 

+ q ^2 Z2P ? ( n i ~ T )dr u 0 + (p+Tj/bii+i.n-F?(«i +1, n) 

T 

- [a + qz 2 + (p + lj/LH.n + c]P*(ni, n) = 0 (19) 

and 


cPt(ri L ,T L - 1) + P^2 z L -iP L (n L - 1, t)S TL t o 

T 

+5'^2P L (.n L -l,T)5 TLt0 + ((3 + q)iJ, nL+1:TL Pl{n L + l,T L ) 

T 

- \pz L —i + 6+ (/3 + q)pn Lt r L + c\Pl(n L ,T L ) = 0. (20) 


In Eq. (18) we recognise the stationarity condition for 
the single site with arrival and departure rates equal to 
(pzi- 1 + qzi + 1 ) and (p +g)/x„, lTl , respectively. Similarly, 
Eq. (19) is the stationarity condition for a single site 
with arrival and departure rates equal, respectively, to 
(a + qz 2 ) and (p + 7)/i ni ,ri, while Eq. (20) has arrival and 
departure rates equal, respectively, toJpzL-i + 5) and 
(/3 + q)p nLjTL . These conditions, in the simplified case 
P"n,r = P’n when r > 0 , allow us to write an approximate 
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stationary distribution for each site l analogous to 
but with modified hopping rates 


Pt(n) = ^U 


, -1 


( 21 ) 


with z c . tl = j2 n zr nr=i w . 


c,i;Z 


and 


Zl 


zl 


a + qz 2 


zi = 



Wc,i-, 1 






Im(c + a + qz2 ) 
c + a + qz 2 + (p + 7 )/m ’ 
pi (c + S +PZ L - 1 ) 
c + 5 +P 2 L -1 + {/3 + q)m ’ 
/ii(c + p 2 i i +qz i+ 1 ) 


c + pzi-i + qzi+i + (p + q)iM ’ 

( 22 ) 

where 1 < l < L. The scenario of Eqs. ( [22| corresponds 
to a one-dimensional lattice where each site l receives a 
uniform particle stream of rate qzi + 1 (pzi_ i) from the 
right (left) neighbour and sends particles according to its 
internal dynamics. The consistency condition of the zi in 
Eq. (22), which is equivalent to the conservation of the 


current along the chain, is satisfied for 


a - 7 Zi = pzi - qzi+i = /3z L - 5 = (j), (23) 


l = 1,2,..., L. The solution of this recursive relation 
yields the fugacity zi and the mean current (j) |31| 


zi 


a/3 (|) [aP-'yS- {a + S){p-q)\ 

/?(§) ('y + p-q)-'y(0-p + q) 

(24) 


«/3 (f) - 7^ 


O') = (p-q ) 


7(P - 9 - /3) + /3(p - g + 7) (|) 


rrr, (25) 


which complete the mean-field solution for the model. 

The approximation results in the separation of the dy¬ 
namics of each site, consequently, the mean-field quan¬ 
tum Hamiltonian can be written as: 


l- 1 

Hrnf = Hieft + ^fright + Hi, ( 26 ) 

1=2 

where H\ 0 f tl .Bright, and Hi are obtained from the generic 
one-site Hamiltonian using the mean-field arrival 
rates. 

While the fugacities ( f24j ) are identical to those of a 
standard ZRP on an open chain [ST] , the effective de¬ 
parture rates w Ctn -i are affected by the time correla¬ 
tions and, significantly, become site dependent. This 
is evident at the level of stationary density and vari¬ 
ance profile, respectively (ni) = zid(\n Z Ct i)/dzi and 
ct/ 2 = (ni 2 ) — (ni) 2 = z 2 d 2 (\nZ c j)/dzi 2 in the mean-field 
approximation. The predicted density profile can be non¬ 
monotonic, which contrasts with the stationary profile of 




TA c = 5 mean-field 

-PA c = 5 mean-field 

- TA c = 0.1 mean-field 

- PA c = 0.1 mean-field 


l 

TA c = 5 simulation 
XXX PA c = 5 simulation 
OOO TA c = 0.1 simulation 
CUtT PA c = 0.1 simulation 


FIG. 4. (Color online) Density profile (ni), and variance 
profile <t; 2 (insets), on a chain of length L = 20 with (a) 
constant and (b) linear departure rates. The results ob¬ 
tained within the mean-field approximation (line vertices) 
are compared with the results computed by means of Monte 
Carlo simulations (markers). TA and PA refer to rates 
(a,/3,j,d,p,q) = (0.2, 0.3, 0, 0,1, 0) and (a, /3,7, S,p, q) = 
(0.1,0.2,0.1, 0.1, 0.55, 0.45) respectively. 


the standard ZRP [31. This feature is indeed present 
in the Monte Carlo simulated density profiles for cer¬ 
tain parameter combinations. In fact, for the parameters 
considered, the agreement between mean-field theory and 
simulation is excellent, except when c is very small (see 
Fig# In Fig. [ 5 ] the mean-field predictions for the per- 
site occupation distributions are checked against simula¬ 
tions. The agreement is again good except for the cases 
with smaller values of c, for which the cross-correlation 
Cij = {(nifij) — (ni)(rij))/{cri<jj) between the occupations 
on site i and j appears to be stronger. This is clear in 
Fig. [6] where we report a negative cross-correlation be¬ 
tween neighbouring-site occupations for small values of c. 
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FIG. 5. (Color online) The mean-field probability distribution (line vertices) of the site occupation numbers in a chain of length 
L = 5 are checked against simulations (markers). Symbols +, x, o, □, A and solid, dotted, dashed, dot-dashed, dot-dot-dashed 
lines of the corresponding color (grayscale) refer to sites l = 1,2,3,4, 5 respectively. Parameter combinations as in Fig. [4] 
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FIG. 6. (Color online) Simulation results for the cross¬ 
correlation Cii,i of on-off ZRP on a chain of length L = 20 
with /j„ = n. (a) PA case. Adjacent sites have negatively- 
correlated occupation numbers. Hopping-rate combinations 
as in Fig. [4] Inset, (b). TA case. Spatial cross-correlations 
appear weaker than the PA case, but with longer range. As c 
grows the correlation is gradually lost and a factorised solu¬ 
tion is realistic. 


In the cases explored above, the values of c have been 
chosen in order to guarantee the existence of a well de¬ 
fined NESS with constant average occupation number. 
In the one-site system seen in Sec. Ill A this choice was 


straightforward, as we can derive exactly the congestion 
threshold. In an extended system with unbounded depar¬ 
ture rates, we expect that any strictly positive value of c 
guarantees the NESS because, although a large number 
of particles can pile up during the OFF phase, they can 
be released arbitrarily quickly during the ON phase. On 
the contrary, the extended system with bounded depar¬ 
ture rates appears to be more interesting. For values of 
c smaller than a certain value, the particles accumulate 
on one or more of the lattice sites. We now compare the 
prediction of the mean-field theory for this congestion 
threshold with the results of Monte Carlo simulations 
performed on a chain of length L = 20. In order to eval¬ 
uate numerically the onset of congestion, we make use of 
the parameter (inspired by (351 ) 


^Hot (^ T At) ^tot(^) f 

At (a + S) 1 


(27) 


where t At and n to t(i) = is the aver¬ 

age total number of particles in the system at time t. 
The parameter k measures the difference between the 
rate at which particles arrive and the rate at which par¬ 
ticles leave the system, scaled with respect to the total 
arrival rate. The congestion occurs when n is strictly pos¬ 
itive [36]. For the one-site model 0 0 with /j,„ = n, 
it is straightforward to show that the expected value 
of an is the positive part of the average growth rate 
a — (3/j,c/(c + a). This allows us to approximate a local 
Ki for the generic site l of a chain, by replacing a and [3 
with the mean-field arrival and departure rates, respec- 
























IV. CURRENT FLUCTUATIONS 
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FIG. 7. (Color online) Congestion transition for the on-off 
ZRP with fi = 1 on a open chain of length L = 20. The 
parameter k obtained by Monte Carlo simulations is plotted 
against c for a TA case and a PA case. Hopping-rate combi¬ 
nations as in Figjdl The mean-field congestion thresholds are 
c mf = 0.4 and ~ 1.4 respectively. These values are pin¬ 
pointed by the mean-field approximated local rez (light lines) 
of the site l where the congestion sets in first. 


tively. The numerical Monte Carlo study of m reveals 
the first site where the congestion sets in, tuning c from 
large to smaller values. In a chain with TA jumps, this 
occurs on site 1, for p < /3, or on the site L, otherwise. 
In the PA case, the congestion can set in on the bulk site 
L — 1, as suggested by the non-monotonic density profile 
of Fig. 0 Not surprisingly, the mean-field theory predicts 
this possibility. 


We define the mean-field congestion threshold c m f as 
the smallest value of c such that none of the sites l of 
the system with Hamiltonian (26) has Ki > 0. In the TA 
case, as long as p > /3, c m f is equivalent to the thresh¬ 
old ci derived in Sec. Ill A for the one-site system with 
boundary rates a and (3. The numerical evaluation of 
k for the whole system, plotted against the mean-field 
estimate in Fig. [7] suggests that c m f is an upper bound 
for the true congestion transition in this case. This re¬ 
lation arises as the TA jumps set the system in a highly 
organised configuration, with wave-like fronts which are 
precursors of the slinky motion observed in Hirschberg 
et al. cun and enhance the particle transport. When 
p < /3, c m f = a 2 /{pp — a) marks exactly the onset of 
the congested phase. Conversely, PA interactions seems 
to promote congestion, as there are many jumps which 
block the site and contribute negatively to the particle 
current. In this case the congestion transition occurs for 
a value of c larger than both Ci and c m f (see Fig. [T|. 


In the next section the out-of-equilibrium aspects of 
this model are further investigated by focusing on the 
fluctuations of the particle currents. 


This section is devoted to the study of the full statis¬ 
tics of the empirical currents ji = Ji/t, where Ji is the 
difference between the number of particle hops from site 
l to site l + 1 and the number of hops from site l + 1 to 
site l. This definition is extended to the input current 
jo and to the output current j l- In order to lighten the 
notation, we simply make use of j and J and explicitly 
specify the bond only when necessary. For t oo, j 
converges to its ensemble average (j). However, for finite 
time it is still possible to observe fluctuations of j from 
this typical value. These fluctuations are quantified by 
means of the scaled cumulant generating function e(s) or 
the rate function e(j). In the following we define these 
concepts. 


A. Large deviation formalism 

Generically, in the long-time limit, the probability 
p(j, t) of observing a current j at time t obeys a large 
deviation principle of the form 

P(j,t) ~ e~ teW . (28) 

To obtain the rate function e(j), we first investigate the 
moment generating function of the total integrated cur¬ 
rent J: 

(e~ sJ ) = (l\e-* J \Pj(t)), (29) 

where J is a diagonal operator whose diagonal elements 
are the set of possible values that the integrated current 
can assume and the probability distribution | Pj{t)) = 
T j P{ n i t j J', t)\ n i T , J) is now defined in the occupa¬ 
tion, clock and current configuration space. The dis¬ 
tribution | Pj(t)) can be obtained from a generic ini¬ 
tial state |P 0 (0)) with J) = (0,0) through | P{t)) = 
e~ Hjt \P o (0 )), where e~ Hjt is the time evolution operator 
in the joint configuration and current space. We diago- 
nalise the operator e~ Hjt by means of a Laplace trans¬ 
form e~ sJ e~ Hjt e sJ . In the configuration subspace, this 
reduces to e~ Ht , where the operator H is obtained mul¬ 
tiplying by e~ s (or e s ) the entries of the original Hamil¬ 
tonian H which produce a unit increase (or decrease) in 
j EZ|. Hereafter, we refer to the tilded operator H as 
the s-modified Hamiltonian. Since (l|e _sJ |Po(0)) = 1, 
then ( e~ sJ ) = (l|e _i ^|P(0)), where |P(0)) now denotes 
a probability vector in the subspace of the occupation 
number and the clock variable. Let us denote by \Pa 0 ) 
the right eigenvector of H associated with the discrete 
smallest eigenvalue A 0 . The long-time limit of the gen¬ 
erating function is accessible through 

{e~ sJ ) ~ (l|P Ao )(P Ao |P(0))e-^‘, t -a oo, (30) 

as long as the pre-factors (1 |Pa 0 ) and (Pa o |P( 0)) are 
finite and a point spectrum exists [see, e.g., [37]]. 
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Although the moment generating function and the con¬ 
jugated variable s have an analogue in equilibrium statis¬ 
tical mechanics, i.e., the Helmholtz free energy and the 
pressure, respectively, they are not as readily accessible 
(we cannot tune s as we can do with the pressure or tem¬ 
perature). However, the generating function (29) helps 
to find out the rate function. In fact, as long the limit 
relation ( |30| ) is valid, we can identify Aq with the scaled 
cumulant generating function (SCGF) 


Equation (34) is hard to solve in general. To gain in- 


e(s) = — lim - ln(e 

t —>OO t 


— sj 


>• 


(31) 


The SCGF in turn gives the convex hull of the rate func¬ 
tion though a Legendre-Fenchel transform [5]: 


e(j) = sup s {e(s) - sj}. 


(32) 


When one of the two pre-factors in Eq. (30) diverges, 


or when the spectrum is entirely continuous, we need to 
employ other methods (see Sec. IVB 3). 


B. Analytical results for the single-site system 

For the single-site system, we study the fluctuations of 
the output current, simply denoted by j. The input cur¬ 
rent can be obtained in the PA case by reflection, while 
in the TA case it is given by a simple Poisson process. 
Despite its simplicity, the single-site ZRP exhibits a rich 
fluctuating behaviour, even in the absence of time cor¬ 
relations [551 - HDj . The introduction of the on-off mech¬ 
anism creates a still more interesting scenario. In fact, 
the study of the fluctuations reveals some aspects of the 
correlations which, in the stationary state, are hidden 
within an effective interaction factor. 


sight into the appropriate structure of p n>s and p 0 N,n,s, 
we study first the simple case with constant departure 
rates. 

Constant departure rates. Let the departure rate be 
p n = p when n > 0. Motivated by the stationary state 
result, we assume here that the factors p 0 N, n ,s and Pn,s 
have no dependence on the occupation number and we 
drop the subscript n with the exception of n = 0 , i.e., 
Pon,o s is distinct from p Q N)S . By direct substitution into 
Eq. (34) we get: 

(38) 

(39) 

(40) 

(41) 


Pa o {n, ON) = p O K,n,spA 0 (n), 
P Ao (n, OFF) = (1 - Pon,n, s )PA 0 (n), 
P Ao {n + 1) = p n+1 , s pA 0 (n). 


(35) 

(36) 

(37) 


-(c + a + <5 - A 0 )(l - Pon.o.s) = 0, 

c(l Pon,0,s) (^ T 3 j 4o)Pon,0,s 

+ (/3e -s + 7)pPon,sPs = o , 
(a + 5e s ) - (c + a + 6 - A 0 )(l - p 0 N, s )Ps = 0, 
c(l - p 0 N,s) - [a + 6 + (/? + 7 )p - A 0 ]pon,s 

+ (/3e~ s + ^)pp 0 N,sPs = 0. 


Equation (38) trivially requires p O N,s ,0 = V while we ex¬ 
pect that p 0 N,s < 1- After a long but straightforward 
algebraic manipulation, the system is solved for 

c + {pe~ s + 7 ) (a + 5e s )/(/3 + 7 ) 


PoN,S — 


c + (/3 + 7 )/r + (/3e s + 7 )(a + 5e s )/(/3 + 7 ) 

(42) 


Ps = 


(a + Se s ) 


(pPon,sY 


An — 


"(1 ~ e<i )- 


(43) 

(44) 


1. Small current fluctuations 

The s-modified Hamiltonian corresponding to the out¬ 
put current is obtained from © multiplying the ladder 
operators and Sa'Y 1 by e~ s and e +s , respectively: 

H = -c(a^ i -g Tl )-a(a)l) i / Tl -ll)-/3(e _s a)( ri d ri -d A r 1 d ri ) 

- 7(«jv! d Tl - d Nl d Tl ) - S(e s aYfJ Tl - 1). (33) 

We concentrate now on the eigenproblem 

(H — A1)\Pa) = 0, (34) 

where \P A ) is the generic right eigenvector and A is 
its eigenvalue. It is convenient to write the eigenvector 
\P Ao ), associated to Aq 7 in a form similar to the station¬ 
ary solution d8l)-([T0|), i.e., with components: 


(P + l) 

^(l-e- s ) + ^ 

/3 + 7 ’ P + 1 

Note that setting s = 0, the factor p 0 N,s becomes the con¬ 
ditional probability P*(ON|n) in the steady state. Also, 
the parameter p s and the eigenvalue A 0 have a coun¬ 
terpart in the stationary probability, in fact for s —> 0 , 
p s -A- zwf 1 , and Ao 0. Consequently, we argue 
that A 0 is the lowest eigenvalue of H and, according to 
Sec. |IV A[ the SCGF at least in the neighbourhood of 
s = 0 . 

For later convenience, we define a modified fugacity 
“ + & ' (45) 


P + 1 

and a modified effective interaction 

w c ,s = PPon,s (46) 

such that p s = z s w~l and w CjS —> p for c — i 00 . It is 
worth noting that, while the bias affects only the fugacity 
in the ordinary ZRP m, it affects both the interaction 
term and the fugacity in the on-off model. 

General departure rates. This paragraph covers also 
the special case with linear departure rates p n = n. Mo¬ 
tivated by the results above, we assume that the com¬ 
ponents of the ground state eigenvector \Pa 0 ) satisfy 
Eqs. (35) (37) with 


Pn,s — z s w c 7 n,si 
^c,n,s — /^nPoN,n,S) 


(47) 

(48) 
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for n > 0. With this assumption, the second row equa¬ 
tion of the eigenproblem (34) is solved for A = A 0 = 
a + /3 — {/3e~ s + 7 )z s and the remaining equations yield a 
solution for z s consistent with (451 and an n-dependent 
effective interaction 

= Mn[c + (/3e~ s + 7 )(a + Se s )/(/3 + 7 )] 

Uc ’ n ’ s (Pe~ s + 7 ) (a + <5e s )/(/3 + 7 ) + c+ (/3 + 7 )/^' 

(49) 

The eigenvalue we obtained is the same as the lowest 
eigenvalue Ho (44) of the s-modified Hamiltonian for the 
standard ZRP |37j . In fact, the affinity between the two 
models appears closer if we work in the reduced state 
space obtained by collapsing the states corresponding to 
r = ON and r = OFF, for each occupation number, and 
considering the sum of their non-conserved probabilities 
Pa 0 (ti) = pA 0 {n, ON) + PA 0 {n, OFF). We notice that 
the vector \P\ ) with components Pa 0 (ti) is the right 


eigenvector with eigenvalue Ho of 
H* = a(a + —l)+S(e s a + —l)+'y(a*-—d*)+/3(e~ 


where 


a. = 


0 

0 

V: 


Wc,l,s 0 

0 Wc,2,s 

0 0 

0 0 


0 

0 

w c , 3,3 

0 


--<), 

(50) 


(51) 




a + = 


/ 0 0 0 ... \ 
10 0 
0 1 0 


V : 


(52) 


■J 


and the operator d* has entries 5ijW Ct i tS ■ The operator 
H* is equivalent to the s-modified Hamiltonian of a stan¬ 
dard ZRP with departure rates w c ^ n ^ s . However, it is not 
a genuine s-modified Hamiltonian for the on-off ZRP as it 
shares only the lowest eigenvalue Ho with H (the higher 
eigenvalues being different, in general) hence it only con¬ 
tains information about the limiting behaviour and does 
not generate the dynamics. 

As a partial conclusion, we underline that both the 
systems with bounded and unbounded rates display the 
fluctuating behaviour seen in the standard ZRP as long 
as the ground state satisfies Eqs. (47) and (48). This is 


certainly true for current fluctuations close to the mean 
(j). However, the effective interaction ui crljS has a de¬ 
pendence on n and s different from the standard ZRP 
and this alters the range of validity of this regime. In the 
following, we show that larger current fluctuations in the 
on-off ZRP can be strongly affected by time correlations. 


2. Range of validity 

The scenario seen so far is an analytical continuation 
of the stationary state. Despite this, certain values of the 


bias s correspond to non-analyticity in the SCGF. Such 
a behaviour is often referred to as a dynamical phase 
transition because of the analogy of the SCGF with the 
Helmholtz free energy. According to Sec. |IV A[ a tran¬ 
sition occurs as soon as the scalar product (1 |Pa 0 ) or 
{Pa 0 |-P(0)) diverges. The choice of the initial distribu¬ 
tion |P(0)) influences the value of the second norm. In or¬ 
der to ensure a finite (Pa„ |P(0)), we will always consider 
an empty site as initial condition, unless explicitly stated 
otherwise. We must also ensure that the norm (Pa 0 \Pa 0 ) 
is finite, i.e., that the eigenvector is normalizable and the 
discrete eigenvalue A 0 exists. We now derive exactly the 
conditions under which the norms (1 |Pa 0 ) and (Pa 0 \Pa 0 ) 
converge and it is possible to identify the SCGF with the 
lowest eigenvalue H 0 given in Eq. (44). 

Linear departure rates. We focus first on the case 
with p n = n. For this particular choice of the inter¬ 
action, particles in the memoryless ZRP can never pile 
up and the current shows a smooth SCGF. On the con¬ 
trary, in the on-off model, the particle blockade alters the 
statistics of small currents. From a mathematical point 
of view, a transition occurs when (l|P/\ 0 ) diverges. The 
condition (l|Pa 0 ) < 00 is satisfied for linij^oo p U:S < 1 
where p UtS is defined in Eq. (47). For later convenience, 
we simplify this condition as 


Ho < c + <5(1 — c s )> 

which is satisfied for s > si, where 

«, 2a/3 


(53) 


a/3 — j35 — (3c — cy + yH a ft 2 5 + {aft — /35 — /3c— cy) 2 

(54) 

In the PA case, si is always finite. In the TA case, i.e., 
( 7 , S ) = ( 0 , 0 ), the critical value si = In [a/{a — c)] is well 
defined only for c < a. 

We can prove that, when p n = n, the norm (Pa 0 \Pa 0 ) 
is always finite. In Appendix [bJ the eigenvector (Pa 0 I 
is derived. Its components have a form similar to 
Eqs. (35)—(371, with the factors p Q NjS and p s replaced 
by Porvi = 1/2 and p\ f ft = (/3e -s + 7 )/(/3 + 7 ) respec¬ 
tively. The series (Pa 0 |Pa 0 ) is simplified by summing 


first the pairs corresponding to the same occupation num¬ 
ber and the condition for convergence can be written as 
lim^^oo p^ps^n < 1, which is always satisfied. Conse¬ 
quently, for linear departure rates, the only mechanism 
responsible for dynamical phase transitions is the on-off 
clockwork, which becomes dominant when (llP^) di¬ 
verges. 

Constant departure rates. Let us consider the case 
p n = 1, n > 0. The scalar product (1 |Pa 0 ) is finite 
when the n-independent parameter p s is less than 1 and 
a dynamical phase transition occurs at p s = 1. In the PA 
case, the solution of this equation for s involves a cumber¬ 
some cubic and therefore is not reported here. However, 
in the TA case, si = ln[a(a — p/3)/(cp/3 — apfi — ca )]. 
In order to check whether (Pa 0 \Pa 0 ) is finite, we again 
need the eigenvector (Pa 0 \- As the dependence on p n 
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cancels in the left eigenproblem, (Pa 0 I is the same as 
the linear departure rate case, see Appendix |Bj The 
condition for convergence is p'f'ps < 1 and the value 
of s such that p'g !t p s = 1 is referred to as S 2 . Also 
here, we only report explicitly the critical bias S 2 = 


— ln[(—c+ \/c 2 + 4c/3)/(2a)] for the TA case. The values 
s = s i and s = S 2 mark the onsets of new phases. 

We notice that the scenario seen so far is entirely en¬ 
coded into the operator H* (50). In fact, this operator 


not only has lowest eigenvalue A 0 , as seen in Sec. IV B 1 


but the normalisation of its ground state eigenvector 
yields sums (l|P^ o ) and (PJ |P£ ) that diverge at the 
same critical points si and S 2 , respectively. In the fol¬ 
lowing, we focus on the large-fluctuation regimes s > s i 
and s < S 2 - 


3. Large current fluctuations 


We employ different approaches to study the large fluc¬ 
tuation regimes in the linear and constant departure rate 
cases. 

Linear departure rates. For this special case, we con¬ 
sider first a finite-capacity version of the on-off ZRP. In 
fact, the SCGF on a discrete finite configuration space is 
always given by the smallest eigenvalue of the s-modified 
Hamiltonian, as the prefactors in (301 are always finite. 


For the TA case, we truncate the Hamiltonian (121 by 


imposing a reflective boundary in the state with occu¬ 
pation number N. The resulting matrix in block form 
is 


H n 


/ C+CK 0 

— C Ot 

0 0 

0 -/3a4i 

0 0 

0 0 


\ 

— OL —ot 

0 0 

C+Q! 0 

— o a+Pni 

0 0 

0 -/3/X2 







0 —fiflN 

V 




c 0 

— e / 


(55) 

which defines a Master equation where the n-th block row 
specifies the dynamics of the configuration with occupa¬ 
tion number n and within each block the first (second) 
row corresponds to an OFF (ON) phase. 

In the present linear departure rate case p n = n and 
the matrix Hn generates the dynamics of a generalised 
exclusion process BU with on-off mechanism, or a queue 
with Markovian arrival times, general service and fi¬ 
nite capacity N [33] [34] • According to the procedure 
of Sec. |IV A[ the finite-capacity s-modified Hamiltonian 
Hn is obtained by multiplying the upper-diagonal rates 
p n ( n = 1,2,3of Hn by e~ s . The numerical 
evaluation of the spectrum of Hn, see Fig. [8j shows that 
the two lowest eigenvalues get closer with increasing val¬ 
ues of N. This gives a clue about the limiting behaviour 



s 


FIG. 8. (Color online) Real part of the spectrum of 
the finite-capacity s-modified Hamiltonian for parameters 
(a, /0,7 ,S,c) = (0.2, 0.3, 0,0, 0.5) and p n = n. For s < si 
the smaller eigenvalue converges to Ao = a(l — e _s ), while 
for s > si it converges to c. 


for TV —> oo, where the eigenvalues coalesce at s = si 
and two different dynamical phases emerge. The SCGF 
converges to a constant branch for s > s±. In the limit 
s — i oo, the truncated s-modified Hamiltonian is lower- 
diagonal and its eigenvalues are given by the escape rates. 
As long as the condition c < a holds, the smallest eigen¬ 
value is c. It corresponds to the escape rate from the 
configuration with N particles and OFF state. We ex¬ 
pect that the corresponding eigenvector does not satisfy 
the ansatz (471 (48). Dynamical phase transitions due 


to the crossover of eigenvectors are observed in spatially- 
extended non-equilibrium models such as the Glauber 
model with open boundaries [32]• We argue that, in the 
infinite capacity limit, the SCGF is given, for s > S \, 
by the escape rate of the system with an instantaneous 
congested state and OFF state. 

Our prediction is checked against numerical simula¬ 
tions, as shown in Fig. 9])a). The simulations employ 
an advanced Monte Carlo algorithm, referred to as the 
“cloning” method, which allows us to measure directly 
the SCGF 1141 143j . This method permits the integration 
of the dynamics generated by an s-modified Hamiltonian 
H 1 by means of the parallel simulation of A f copies of the 
system. A system in state i may be cloned or pruned 
with exponential rate Hu — y~h H,j , in order to account 

for the fact that H does not conserve the total proba¬ 
bility. The average cloning factor gives the SCGF. This 
prescription is believed to be exact for Af —> oo, t —> oo, 
and is not reliable when the cloning factor is larger than 
J\f (shaded areas in Fig. [9] and 121, as studied in \M\ . 
Our implementation correctly reproduces the most rele¬ 
vant features of the SCGF, i.e., the non-analyticity in si 
and the constant branch for s > s i, but loses accuracy 
for large positive currents (s < 0) presumably due to the 
finite J\f effect. 
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FIG. 9. (Color online) SCGF of the on-off ZRP 

with fj, n = n for (a) TA hopping rates, (ct, /3, 7 , <5, c) = 
(0.2,0.3,0, 0, 0.1), and (b) PA hopping rates, (a,/3, 7 , <5, c) = 
(0.1, 0.2,0.1, 0.1, 0.1). Points are data from the cloning al¬ 
gorithm, A f = 10 4 ,t = 10 4 . The SCGF is systematically 
overestimated for small values of s. We expect a better ap¬ 
proximation but a slow convergence for larger ensembles and 
longer simulation times. 


In the PA process, the lowest eigenvalue does not ap¬ 
pear to converge to a finite value in the limit s —> 00 . 
From the condition (53) for the eigenvalue crossover, we 
suggest 


e(s) 


^(l-e- s ) + ^(l-e s ), S < Sl 
c + <5(1 — e s ), s > si. 


(56) 


The right branch can be physically understood by sepa¬ 
rating the contributions of the particles leaving the site 
rightwards, which contribute a term c as in the TA case, 
and the particles injected from the right boundary, which 
independently follow a Poisson process with rate 5 and 
contribute a term d(l — e s ). Since in this regime the 
particles pile up, the corresponding SCGF branch does 
not depend on the left boundary. Numerical simulations, 
shown in Fig. |9])b), confirm our argument. There is no 
analogue, in the menrory-less ZRP, of the c-dependent 


dynamical phase for s > s 1 , which arises as a conse¬ 
quence of the temporal correlations. 

For SCGFs with non differentiable points, as in 


Eq. (56), the Legendre-Fenchel transform (32) of e(s) 
gives m general the convex hull of the rate function e(j), 
which can hide a non-convex shape. However, for this 
system, we argue on physical grounds (see following) that 
Eq. (32) gives indeed the true rate function, i.e., 


e(j) 


= 


z + S + j 

~ j ln(- 

?). 

j < 

jl,a 

-Slj + c 

+ <5(1 - 

-e fll ), 

jl,a 

< j 

a/3 

7<5 

Ji 2 | i 'i 5 



/3 + 7 

/3 + 7 

V /3 + 7 /3 + 7 




j + 1 

/,'2 1 4 aH I s 

J /3+7 /3+7 

j> 

31 ,b- 

+j 

In- 

2 #- 




J /3+7 

(57) 

The two critical currents j i iQ = — 8e Sl and j-\ 7 = 
j^e~ Sl — are, respectively, the right and left 

derivatives of e(s) at s = S\. In the TA process ji^ a = 0. 
The phase j < j i j0 is obtained from the Legendre-Fenchel 
transform of e(s) in the interval s > s 1 , while the phase 
j > j 1,6 is derived from e(s), with s < s\. The transition 
value Si is mapped to the linear branch in j-[ a < j < 

This behaviour is equivalent to an ordinary equilibrium 
first-order phase transition, where a linear branch of a 
thermodynamic potential corresponds to the coexistence 
of two phases. In this non-equilibrium system, the mixed 
phase consists in a regime where, for some finite frac¬ 
tion of time, the current assumes value j i jQ , while for the 
rest of the time it has value ji .5. As a result, the rate 
function in this region is linear with j, as predicted by 
the Legendre-Fenchel transform. This argument is sup¬ 
ported by standard Monte Carlo simulations (ensemble 
size of 10 10 ), and it is particularly evident in the TA case 


(Fig. [10|. 

The different phases can be physically understood by 
observing the effect of the particle blockade. In the case 
with TA hopping rates, when the site is OFF, the parti¬ 
cles accumulate and the outgoing current is necessarily 
zero. The zero current is mapped to the flat section of 
the SCGF. This is the dominant mechanism responsible 
for zero current. At the end of an OFF period, we have a 
configuration with many particles on the site. When the 
lock is released, particles can leave the site with a rate 
proportional to the occupation number. Consequently, 
the particles are quickly released after an OFF period 
and the current jumps to a positive value. In particu¬ 
lar, the probability of having currents larger than j\ j, 
is dominated by the phases in which the site is ON. In 
the presence of arrivals from the right boundary (d 7 ^ 0 ), 
the blocked configuration becomes important for negative 
currents j < ji, Q , and the rate function has an additional 
term corresponding to an independent Poisson process 
with rate <5. 

As an aside, the dynamical phase transition seen at Si 
is not restricted to the particular on-off ZRP explored 
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the PA counterparts of the effective departure rate and 
fugacity found in Sec. Ill A while the limit x —> 0 cor¬ 
responds to the empty-site state. The exact calculation 
of the full spectrum and of its eigenvectors, reported in 
Appendix [Cj gives the following representation: 


<l|e-^|P(0)> 


i-x r e-^n 
27 TiX(j)J Cx (C- ^)(C- l) 


1 — x f (y£ — l)e 

27 rix fc 2 (c - ^)(C - Me - y) 


(58) 


where e(£) is obtained from the expression for the con¬ 
tinuous band of the spectrum e(k) after the substitution 
C = e ik and 


FIG. 10. (Color online) Rate function e(j) for the on-off TA 
process with (a,/J, 7 , <5, c) = (0.2,0.3,0,0,0.01) and = n 
(solid line). The points are simulation data for the finite¬ 
time rate function — ln[Prob(J/t = j)\/t computed at t = 
200, 300,400, 500, 600, 700 (top to bottom). 


here. For example, an alternative on-off ZRP with un¬ 
bounded departure rates and on-off dynamics indepen¬ 
dent from the arrivals, displays the same fluctuating sce¬ 
nario. Also, spatially extended spin systems such as 
the contact process [14] and some kinetically constrained 
models ]45| , can possess active and inactive phases coex¬ 
isting at s = 0. 

Constant departure rates. In this case when fi n = 1, 
n > 0, the operator H has a continuous band which 
governs the fluctuations in certain regimes. A way to 
obtain the SCGF is to evaluate the long-time limit of 
the matrix element (l|e _ ' ff *|P(0)) by computing the full 
spectrum and the complete set of eigenvectors of H. This 
task appears to be rather complicated for the s-modified 


(/3e s + 7)w e ,s 
(a + <5e s ) ’ 


(59) 


Hamiltonian (33), requiring spectral theory and integral 
representation of block non-stochastic operators [46]. As 


an approximation, we can use the reduced operator (50) 
and study the simpler expectation {\\e~ H *|P(0)). Recall 
that H* has the same lowest eigenvalue Aq as H, at least 
in the regime S 2 < s < Si where the ansatz (47) (48) 


is valid. Outside this regime it is expected to yield only 
approximate information about the current fluctuations. 

The integral representation allows us to take into ac¬ 
count the dependence of the fluctuations on the initial 
condition. We follow the same procedure as [38] [39], 
with the difference that the departure rate here depends 
on s. In fact, the solution found only has a weak depen¬ 
dence on the functional form of ru CjS but, nevertheless, we 
report the explicit calculations for completeness. As ini¬ 
tial condition, we choose a geometric distribution with 
parameter x, i.e., |P(0)) = (1 — %) xn \ n ) where 
| n) denotes the configuration of the site with n parti¬ 
cles and is an element of the natural basis for H*. The 
steady state is obtained for x = zw ~ l , where w c = 
n(a+6+c)/(a+5+c+((3+'y)iJ,) and z = (a+6)/U3+ , y) are 


V = (o , \ - J(a + Se s )(Pe s + 7 )w Cj8 , (60) 

(p + l)w c ,s v 

e(k) =a + S + (/3 + 7 )ui C jS 

— 2y/(a + 5e s )(/3e~ s + 7 )w c ,s cos(fc). (61) 

The integration contours C\ and C 2 are anti-clockwise 
circles centred around the origin with radius tfi^ 1 < |£| < 
{(j)x)- 1 and infinitesimal size respectively. 

The long-time limit of this integral is computed by 
means of the method of steepest descents with saddle 
point at C = 1- When the saddle point contour engulfs 
one of the poles of the integrand we must also take into 
account the residue m- 

For fixed parameters s and x, the leading term in the 
long-time limit of (l|e _ ' ff **|P(0)) is given by the slowest 
decaying exponential and the SCGF is determined by one 
of the rates e{y), e(l), e((x^) _1 ). Tuning s or x, the 
positions of the poles with respect to the saddle point 
contour are altered and the leading term in the integral 
expansion changes. This produces the phase diagram of 


Fig. 11 for the SCGF. The critical line corresponding 
to the solution of </> = y is s = Si. The line s = S2 
corresponds to y = 1. These two phase transitions were 
also found in Sec. IV B 2| as critical points in the full state 
space. The curves s = S3 and s = S4 are solutions of 
1 = (</>£)” 1 and y = (<j>x)^ 1 respectively. The tri-critical 
point S 3 = S 4 is at x c . It is worth noting that higher 
positive current fluctuations retain a dependence on the 
initial condition x and that, unlike the memoryless ZRP, 
the critical point si can fall in the positive current range. 
The explicit expressions in terms of s for the TA case are 
reported in Appendix |D| 

We distinguish four phases: 

Phase A: s > si. In this case the leading term arises 
from the pole at £ = (/>. The product (1|Pa 0 ) diverges 
and the SCGF is different from the lowest eigenvalue Aq, 
being given instead by 

e(s) = <5(1 - e s ) + /3w CtS ( 1 - e~ s ). 


(62) 
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FIG. 11. (Color online) Phase diagram, based on Eq. ( |58| ), 
for the current fluctuations in the PA process with fi n = 1 
and (a, /3, 7 , 5, c) = (0.1, 0.2,0.1, 0.1.0.5). The lines Si and S 4 
correspond to first-order dynamical phase transitions while S 2 
and S 3 mark second-order transitions. 


This phase corresponds to very small positive currents 
(in particular when S = 0) or large backward currents. 
Large negative currents are mainly governed by the rate 
S of particle arrival from the right, which contributes to 
the SCGF with the first term of (62). The second term 


corresponds to particles that jump rightwards from the 
site with an effective rate /3w C}S . The current fluctua¬ 
tions in this phase are optimally realised by a site with 
arbitrarily large occupation number (instantaneous con¬ 
densation) that acts as a reservoir, so that the outgoing 
current has no dependence on the left boundary hops [39]. 
We argue that the presence of a left and a right term in 
Eq. (62) is generic for this phase, although there is no a 


priori reason for the effective rate w C}S to have the same 
form as in the small fluctuation regime. In the PA case, 
for large values of s, the SCGF is dominated by the first 
term and is not sensitive to the functional form of w r 


Phase B: ((s 2 < s < Si) A (x < x c )) V ((S 4 < s < 
si) A (x > x c ))). This phase arises when the pole at 
C = y, corresponding to the lowest eigenvalue A 0 (441, 
becomes dominant, hence 


e(s) = - 5 —(1 — e _s ) + -p—— (1 — e s ). (63) 


/3 + 7 


/3 + 7 


The probability of fluctuations in this regime is asymp¬ 
totically identical to the standard ZRP. In this range 
the site has finite occupation and the probability that 
a particle leaves is conditioned to an arrival event, just 
as in 


Phase C: {x < x c ) A (S 3 < s < s 2 ). This phase arises 
from the saddle-point at £ = 1. It corresponds to a 
large forward current sustained by a large inward cur¬ 
rent from the left boundary. The asymptotic form (30) 


FIG. 12. (Color online) SCGF of the on-off ZRP with 
(a, /?, 7 , 5, c) = (0.1, 0.2, 0.1, 0.1, 0.5) and p n = 1. Points are 
data from the cloning simulations, A f = 10 4 ,t = 10 4 . Dotted 
line is the SCGF of the Markovian-ZRP (c —► 00 ) with same 
boundary rates. Solid line is the analytic approximation (621- 


(64). The SCGF of the ZRP with s-independent departure 
rate w c would overlap the solid line at this scale. 


ground state. This also represents an instantaneous con¬ 
densate, but with particle number growing as the square 
root of time [33. Here the spectrum of H * is continuous 


and the SCGF is given by the minimum of the band (61): 


e(s) = a + S + (/3 + 7 )w c>s - (a + 5e s )(/3e s +7)iu c , s . 

(64) 

Phase D: [(s < S 3 ) A ( x < x c )\ V [(s < S 4 ) A (x > x c )\. 
This phase arises when the residue at £ = (ifix)- 1 domi¬ 
nates the long-time behaviour: 

e(s) = a+S+(/3+'y)'W CiS -(Pe~ s +"f)w CjS x-(a+Se s )/x. 

(65) 

It corresponds to a large forward current of particles that 
is most likely to be realized from an initial configuration 
with very high occupation number and also has an ana¬ 
logue in the standard ZRP [59] , 

These results are compared to the cloning simulations 
in Fig. [12] for x —> 0. Similarly to the independent- 
particle case, the cloning data for the left branch, cor¬ 
responding to large positive currents, is potentially af¬ 
fected by finite- AT effects [33]. It turns out that for the 


chosen parameters our approximation (62)—(64), plotted 


as a solid line, is very close to the naive approach (not 
shown) in which the same representation ( |58[ ) is used, but 
the effective departure rate has the s-independent form 


w c (see Sec. Ill A) for all the regimes. The analytical 


SCGF does not match the simulation points in either of 
the phases A and C. We attribute this to the failure of 
the assumption (47) (48) for the ground state in phases 


still holds, but with an oscillating (non-decaying in n) 


A and C. In other words, large fluctuations cannot be ex¬ 
actly described by an effective departure rate w CtS with 
a simple functional dependence on s. 
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FIG. 13. (Color online) Rate function for the on-off ZRP 
with (a, /3, 7 , 5, c) = (0.1,0, 2, 0.1, 0.1, 0.5) and p n = 1. Points 
are data for — ln[Prob( J/t = j)\/t from standard Monte Carlo 
simulation at times t = 100, 200,300, 400, 500,1000, 2000. 
(top to bottom). The solid line is the analytical approxi¬ 
mation for t —^ oo. 


FIG. 14. (Color online) Simulation results for the SCGF in 
a five-site on-off ZRP with /i n = 1 and (a, /3,'y,5,p,q,c) = 
(0.1,0.2,0.1,0.1,0.55,0.45,0.5). The solid line is the c- 
independent expression for the lowest eigenvalue of the s- 
modified Hamiltonian for the five-site Markovian-ZRP l39| . 


In Fig. 13 the rate function e(j), computed by means 
of a Legendre-Fenchel transform on the SCGF (62) (64), 
is compared to the finite-time rate function obtained from 
standard Monte Carlo simulations with an ensemble size 
of 10 10 . Although approximate, e(j) appears to capture 
well the shape of the long time limit for the simulation 
data points. 


C. Numerical results for large system 

The lack of a stationary product form solution for the 
on-off ZRP on an extended lattice makes the analytical 
study of fluctuations, across the generic bond, impracti¬ 
cal. It would be possible to use the mean-field stationary 
solution to derive an approximate SCGF using the same 
procedure as in the single-site model. However, we do 
not expect the result to be accurate, especially for small 
values of c and for current fluctuations far from the mean. 
To explore the larger system we make use of the cloning 
method, see Fig. [M] While the statistics of rare currents 
is bond-dependent, it is possible to appreciate that for 
each bond the SCGF matches that of a Markovian ZRP 
in the neighbourhood of s = 0 , a feature shared with the 
one-site system. 

The central regime satisfies a Gallavotti-Cohen fluc¬ 
tuation symmetry [H[] e(s) = e(E — s ), with E = 
ln[(p/< 7 ) i- 1 a/ 3 /( 7 < 5 )]. Such a relation seems to be en¬ 
sured by the fact that the relative probabilities of parti¬ 
cle jumps towards the left or the right are independent of 
the time that the particle spends on a site. This property 
is related to the direction-time independence of Ref. |50| . 
However, the fluctuation symmetry is not guaranteed to 
hold on an arbitrary domain in systems with infinite state 
space [55] , In fact, as expected, we see here a c-dependent 


breakdown for large fluctuations. 


V. DISCUSSION 

We have studied an open-boundary zero-range pro¬ 
cess that incorporates memory by means of an additional 
“phase” variable. The particles are blocked on a lat¬ 
tice site (“phase OFF”) when a new particle arrives and 
consequently congestion is facilitated. After an expo¬ 
nentially distributed waiting time with parameter c, the 
block is removed (“phase ON”). At first sight, the effects 
of time correlations are hidden. The stationary state so¬ 
lution of the one-site system can be written as in the 
Markovian case, with an effective on-site interaction w Cjn . 
This means that, if the direct interactions are unknown, 
it is not possible to distinguish a single site with on-off 
dynamics from a standard memoryless ZRP by looking 
only at the occupation distribution. 

However, the presence of ON and OFF phases alters 
the statistics of the outwards particle hops. This be¬ 
comes important in the spatially extended system where 
each site receives particles, from its neighbours, accord¬ 
ing to a non-Markovian process. As a consequence, a 
product form solution is in general not expected and we 
have relied on a mean-field approach for the analytical 
treatment. This approximation consists of replacing the 
true particle arrival on each site with a memoryless pro¬ 
cess, while keeping exact information about the on-site 
particle departure as well as the lattice topology. This 
procedure can be applied in principle to decouple non- 
Markovian ZRPs on an arbitrary lattice, provided that 
it is possible to solve the consistency equation for the fu- 
gacities. We found that, in the chain topology studied 
here, the mean-field approach is very accurate for large 
values of c and gives an analytical estimate c m f for the 
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congestion threshold. 
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tions close to the mean current are obtained by analytic 
continuation of the stationary state and are indistinguish¬ 
able from the fluctuations in a memoryless ZRP. How¬ 
ever, under certain conditions large current fluctuations 
are optimally realized by the instantaneous piling up of 
particles on the site and the statistics of such fluctuations 
change abruptly. In the absence of direct inter-particle 
interaction we have found a memory-induced dynamical 
first-order phase transition, i.e., the scaled cumulant gen¬ 
erating function (SCGF) e(s) is non-analytic at a partic¬ 
ular value s'i. In the totally asymmetric case, this occurs 
only if the parameter c is smaller than the arrival rate a. 
The system with constant departure rates, i.e., attractive 
inter-particle interaction, undergoes second-order as well 
as first-order dynamical phase transitions. The state of 
the system during a small fluctuation event has the same 
form as the stationary state, but with a more general 
modified effective interaction factor. Indeed, the exact 
phase boundaries and the large deviation function of this 
regime are encoded in the reduced operator H * [Eq. (50)], 
which has the same structure as the s-modified Hamil¬ 
tonian of the standard ZRP, but with an s-dependent 
effective interaction factor. We have used the same op¬ 
erator H* to find an approximate solution for the fluc¬ 
tuations outside this phase. Numerical tests confirm the 
presence of the predicted c-dependent dynamical phase 
transitions. 


The separation between a small-fluctuation regime, 
with a memory independent SCGF, and high-fluctuation 
regimes, where memory plays a more obvious role, is a 
feature also found numerically in the spatially-extended 
system. It would be of interest to explore the role of 
topology in more detail as well as to look for similar 
memory effects in other driven interacting-particle sys¬ 
tems. Furthermore, we point out the importance of solv¬ 
ing the eigenproblem (34) for the full s-modified Hamil¬ 
tonian H [Eq. 


which provides exact information 
about the strongly fluctuating regimes. This would be 
of interest in queueing theory; in fact quasi-birth-death 
processes, which contain as a special case the single-site 
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Appendix A: Derivation of the stationary state 

Summing Eqs. Q and ©>, and imposing the station- 
arity condition, it follows that 

/3/z„ + iP*(n + 1, ON) - aP*(n ) 

= /3p, n P*(n, ON) — aP*{n — 1), (Al) 

while the stationarity conditions on Eqs. <[5j) and ([7| im¬ 
ply the boundary conditions 

/?/xiP*(l, ON) — aP*(0, ON) = 0, (A2) 

P(0, OFF) = 0, (A3) 

which, together with ( |Al| ), allow us to write the recursive 
relation 

Pii n+l P*(n+ l,ON) = aP*(n, ON) + aP*{n, OFF). 

(A4) 


Using the stationarity condition on Eq. © 

(a + c)P*{n + 1, OFF) = aP*(n, ON) + <xP*(n, OFF), 

(A5) 


we eliminate P{n) from Eqs. (A4) and (A5) and get 


(a + c)P*{n + 1, OFF) = /3/j, n+1 P*(n + 1, ON), (A6) 


hence, 


P*(n, OFF) 
P*(n, ON) 


a + c + Pfi 
(a + c ) 
a + c + /3 /j, 


~P*(n ), (A7) 

n 

— P* (n). (A8) 

n 


The ratios (/I/x n )/(a + c + n n /3) and (a + c)/(a + c + 
Pun) are the conditional probabilities P*( OFF n) and 
P*(ON|n), respectively. Substituting in (A4) or (A5) we 
get the recursive relation: 


M“ + , c) p. ( „ + 1) 

a + C + PHn 


a 

1 


P*(n). 


Finally, iterating and using the definitions of z and Z c 
we find the probability mass Q. 
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Appendix B: Left eigenvectors of the s-modifled 
Hamiltonian 

The derivation of (Pa 0 I when y n = p, n > 0, is as 
follows. Assuming that the left-eigenvector components 
satisfy 


P lGft (n, ON) = p^ iS P left (n), (Bl) 

Pie ft (n, OFF) = (1 - JPieftM, (B2) 
fleft(n+l) =pTPieit(n), (B3) 

we get the explicit equations 

— (a — A + c +$){1— Pon,s,o) + 

+ (a + 6e s y s eft ( 1 - Po f N,J = 0, 

(B4) 

(a + Se s ) P T (1 - pgj, J - (a - A + s , 0 = 0, 

(B5) 

- (a - A + c + <5)(1 - Pon,J + c (1 - Pon, s ) 

+ (a + <5e s ) Pg 6ft (l - PoS.J = 0. 

(B6) 


(/3e s + 7)^“, s ,o - (a - A + (/3 + y)/z + 

+ (a + 5e s ) p l f 2 [l - Pon, s ) = 0, 

(B7) 

-pfVot.Ja - A+p(fi + 7 ) + 5} + np'°* tS (7 + /3e~ s ) 

+ P r 2 (i-PoN lS )( a +^ s )=o, 

(B8) 


where the factor s 0 is assumed to be different from 
-Pon s by analogy with the right eigenproblem. The 
Eqs. ( |B4j ) and (B5) give (a-A + c + (5)(l—Po“ jSi0 )-|-(a — 
A + &W 0 S, s ,o + cpAn,6.o = °> which is verihed forp^“ s0 = 
§• The Eqs. (B7} and (|B8]) imply p^ s = p‘'"; 6 . 0 . Af¬ 
ter the substitution, the remaining equations are solved 
for A = A 0 and p 1 ®* 4 = (/3e _s + 7)/(/3 + 7 )- With t hose 
constants, it is easy to verify that the ansatz (Bl) (B3) 
is consistent even in the general departure rate case. In 
fact, after substitution, all the terms containing p n can¬ 
cel out. In the reduced state space we get a consistent 
result since the row vector (PJ | with components given 

by ® satisfies (PT \H* = A 0 (PT|. 


Appendix C: Spectrum and integral representation 


where (m| (|n)) is a row (column) vector with a “1” in the 
m-th (n-th) position and “0” elsewhere. To evaluate the 
right-hand side of ( |C1| ), we first seek for normal modes 
of the dynamics generated by the operator H* ( |50| ). We 
transform H* into the symmetric form , where 

$ is the diagonal operator with entries Sijtf) 1 , Sij is the 
Kronecker delta, i,j = 0,1, 2,..., and <fi is the combina¬ 
tion of parameters (59) in the main text. The associated 


eigenproblem is solved after a Fourier transformation. Its 
eigenvalue e(fc) [Eq. (61)], has eigenvector \ij)'{k)) with 
components 2/7T sin(n/c -+• <p). Substituting this in the 
first row equation for the eigenproblem, we get the fol¬ 
lowing expression for <p: 


e i2 V = 


1 - e^y 


1 — e 


—ik 


y 


(C2) 


where y is given by the s-dependent expression (60). For 
y < la discrete eigenvalue appears with eigenvector 
\tp'{0)) = \J\ - y 2 E)T=o V n \ n ) and eigenvalue A 0 [Eq. 
(44)] while, for y > 1, the infimum of the spectrum is 
given by e(0). 

The vectors \^{k)) = ^~ 1 \ip'(k)} , k £ (0,7r], 

and |^(0)) = $ _1 |^/> / (0)) form a complete set, i.e. 
fg | i/j(k))(ip(k)\dk+ |^(0))(V’(O)I = 1. Inserting this rep¬ 
resentation of the identity in Eq. (Cl), the right-hand 
side becomes 

CO /--7T 

(1 — x) ^ x n / (m\il>(k)}('i/j(k)\n)e~ e( - k ^ t dk 

n,m =0 ® 

co 

+ 0(1 — y)(i — x) x n r- m (i-y 2 )y n+m e- Aot , 

m,n=0 

(C3) 

where 0 denotes the Heaviside step function. Using the 
fact the eigenvectors are odd in k , the integral in Eq. (C3) 
can be rewritten as 

J ^ e ik(n-m) _ e ik(n+m) e i2^ e ~e(k)t d ^ ( C4 ) 


and, using Eq. ( |C2[ ), it becomes 

c 


ICI=i 


n—m—1 _ ^-n-\-m ^ ClJ \ 


C -y 


'j e- £(C)t dC, (C5) 


where £ = e lk and e(C) = e[/c(<C)]. Deforming the inte¬ 
gration contour to Ci for the first term in the integrand 
and to C 2 for the second term, we obtain the representa¬ 
tion (58). The last term in Eq. (C3) cancels out with a 


pole contribution at £ = y for y < 1. 


In this appendix, we report the calculations which lead 
to the integral representation (58). Let us impose an 
initial condition of Boltzmann type for the system, so 
that 


Appendix D: Phase diagram for the current 
fluctuations in the TA process with bounded 
departure rate 


OO 

<l|e-^|P(0)) = (l —z) £ x n (m\e~ nH \n), (Cl) 

n,m =0 


In this appendix, we report the analytical forms of 
the c-dependent transition lines between the dynamical 
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phases of the TA case with p n = n. The resulting phase 
diagram is similar to the PA case (Fig. [TTj), but with the 
transition line identified by Si mapped to a positive value 
of the current. 


s — si: The knowledge of \Pa 0 ) is sufficient to verify 
when the pre-factor (s|P 4 0 ) is finite, i.e., 


a(c + 
(c + 

s < Si = In 


e S a + Pn) 
e~ s a)/3ii 

( a(pn -a) \ 

\ ca — c(3/j, + a/3p J 


(Dl) 

(D2) 


Notice that this condition makes sense when the denomi¬ 
nator in the argument of the logarithm in (D 2 ) is positive, 
i.e., c < a/3p/(ffp — a), while the stationarity condition 
a < /3w c ensures that the numerator is positive. The 
phase boundary can also be obtained from solving <f> = y. 


the region where the condition p'f'ps < 1 holds, i.e., 


a ( c + e s a + f3p) c _ s 1 

(c + e~ s a)f3p 


s > S 2 = — In 


/ \/c 2 + 4c/i 
^ 2 a 



(D3) 

(D4) 


It corresponds to a solution of y = 1. 

s = s 3 : This line corresponds to (cfx)^ 1 = 1. The criti¬ 
cal point S 3 satisfies 


0 -«3 _ 


\J (cpx 2 + a 2 ) 2 + 4,a 2 p 2 x 2 — cpx 2 + < 
2apx 2 


(D5) 


s — s 4 : This phase boundary is c-independent, specifi¬ 
cally 

e " 84 = a; -1 . (D 6 ) 


It corresponds to the condition y = ( <f>x) 1 

s — S 2 : This critical point marks the left boundary of 
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